Hands-on-Exercise 3.1: 1st Order Spatial Point Patterns Analysis Methods

Published

January 17, 2024

4.1 Overview

Spatial Point Pattern Analysis is the evaluation of the pattern or distribution, of a set of points on a surface. The point can be location of:

  • events such as crime, traffic accident and disease onset, or
  • business services (coffee and fastfood outlets) or facilities such as childcare and eldercare.

Using appropriate functions of spatstat, this hands-on exercise aims to discover the spatial point processes of childecare centres in Singapore.

The specific questions we would like to answer are as follows:

  • are the childcare centres in Singapore randomly distributed throughout the country?
  • if the answer is not, then the next logical question is where are the locations with higher concentration of childcare centres?

4.2 The data

To provide answers to the questions above, three data sets will be used. They are:

  • CHILDCARE, a point feature data providing both location and attribute information of childcare centres. It was downloaded from Data.gov.sg and is in geojson format.

  • MP14_SUBZONE_WEB_PL, a polygon feature data providing information of URA 2014 Master Plan Planning Subzone boundary data. It is in ESRI shapefile format. This data set was also downloaded from Data.gov.sg.

  • CostalOutline, a polygon feature data showing the national boundary of Singapore. It is provided by SLA and is in ESRI shapefile format.

4.3 Installing and Loading the R packages

In this hands-on exercise, five R packages will be used, they are:

  • sf, a relatively new R package specially designed to import, manage and process vector-based geospatial data in R.

  • spatstat, which has a wide range of useful functions for point pattern analysis. In this hands-on exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.

  • raster which reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this hands-on exercise, it will be used to convert image output generate by spatstat into raster format.

  • maptools which provides a set of tools for manipulating geographic data. In this hands-on exercise, we mainly use it to convert Spatial objects into ppp format of spatstat.

  • tmap which provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.

Use the code chunk below to install and launch the five R packages.

install.packages("maptools", repos = "https://packagemanager.posit.co/cran/2023-10-13")
package 'maptools' successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\geral\AppData\Local\Temp\RtmpiGSn1E\downloaded_packages
pacman::p_load(sf, raster, spatstat, tmap)

4.4 Spatial Data Wrangling

4.4.1 Importing the spatial data

In this section, st_read() of sf package will be used to import these three geospatial data sets into R. The 3 data are childcare (preschoollocation), coastal outline, and MP14 subzone.

childcare_sf <- st_read("data/childcareservices.geojson") %>%
  st_transform(crs = 3414)
Reading layer `ChildCareServices' from data source 
  `C:\glimjw\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data\ChildCareServices.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
sg_sf <- st_read("data/MasterPlan2019SubzoneBoundaryNoSeaGEOJSON.geojson") %>%
  st_transform(crs = 3414)
Reading layer `MasterPlan2019SubzoneBoundaryNoSeaGEOJSON' from data source 
  `C:\glimjw\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data\MasterPlan2019SubzoneBoundaryNoSeaGEOJSON.geojson' 
  using driver `GeoJSON'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY, XYZ
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
mpsz_sf <- st_read(dsn = "data/geospatial", 
                layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\glimjw\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

Before we can use these data for analysis, it is important for us to ensure that they are projected in same projection system.

DIY: Using the appropriate sf function you learned in Hands-on Exercise 2, retrieve the referencing system information of these geospatial data.

st_crs(childcare_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(sg_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(mpsz_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

DIY: Using the method you learned in Lesson 2, assign the correct crs to mpsz_sf and sg_sf simple feature data frames.

The cr information isn’t appropriate. childcare_sf  and ‘sg_sf’ are in WGS84, while the ‘other two are’mpsz_sf’ is in SVY21

mpsz_sf <- st_transform(mpsz_sf, crs= 3414)
st_crs(mpsz_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

4.4.2 Mapping the geospatial data sets

Let’s plot a map to show their spatial patterns for each dataframe

childcare_sf consists of spatial points, so it cannot accept tm_fill/tm_borders/tm_polygons. So we have to use tm_dots()

invalid_mpsz_sf <- which(!st_is_valid(mpsz_sf))
invalid_sg_sf <- which(!st_is_valid(sg_sf))
invalid_childcare_sf <- which(!st_is_valid(childcare_sf))

# Print the indices of the invalid geometries
print(invalid_mpsz_sf)
[1]  19  20  24 122 123 128 258 302 320
print(invalid_sg_sf)
[1] 35 59
print(invalid_childcare_sf)
integer(0)
# Now, try plotting again
tm_shape(sg_sf) +
  tm_polygons() +
tm_shape(mpsz_sf) +
  tm_polygons() +
tm_shape(childcare_sf) +
  tm_dots()

Alternatively, we can prepare a pin map by using this code chunk:

tmap_mode('view')
tm_shape(childcare_sf)+
  tm_dots()
tmap_mode('plot')

Notice that at the interactive mode, tmap is using leaflet for R API. The advantage of this interactive pin map is it allows us to navigate and zoom around the map freely.

4.5 Geospatial Data wrangling

Although simple feature data frame is gaining popularity again sp’s Spatial* classes, there are, however, many geospatial analysis packages require the input geospatial data in sp’s Spatial* classes. In this section, you will learn how to convert simple feature data frame to sp’s Spatial* class.

childcare <- as_Spatial(childcare_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(st_zm(sg_sf))

Show the information of these 3 Spatial* classes

childcare
class       : SpatialPointsDataFrame 
features    : 1925 
extent      : 11810.03, 45404.24, 25596.33, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :    Name,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    Description 
min values  :   kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>100044</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>44, TELOK BLANGAH DRIVE, #01 - 19/51, SINGAPORE 100044</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td></td> </tr><tr bgcolor=""> <th>NAME</th> <td>PCF SPARKLETOTS PRESCHOOL @ TELOK BLANGAH BLK 44 (CC)</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>349C54F201805938</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20211201093837</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
max values  : kml_999,                                            <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>99982</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>35, ALLANBROOKE ROAD, SINGAPORE 099982</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td></td> </tr><tr bgcolor=""> <th>NAME</th> <td>ISLANDER PRE-SCHOOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>4F63ACF93EFABE7F</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20211201093837</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
mpsz
class       : SpatialPolygonsDataFrame 
features    : 323 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 15
names       : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C,       REGION_N, REGION_C,          INC_CRC, FMEL_UPD_D,     X_ADDR,     Y_ADDR,    SHAPE_Leng,    SHAPE_Area 
min values  :        1,          1, ADMIRALTY,    AMSZ01,      N, ANG MO KIO,         AM, CENTRAL REGION,       CR, 00F5E30B5C9B7AD8,      16409,  5092.8949,  19579.069, 871.554887798, 39437.9352703 
max values  :      323,         17,    YUNNAN,    YSSZ09,      Y,     YISHUN,         YS,    WEST REGION,       WR, FFCCF172717C2EAF,      16409, 50424.7923, 49552.7904, 68083.9364708,  69748298.792 
sg
class       : SpatialPolygonsDataFrame 
features    : 332 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :   Name,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       Description 
min values  :  kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>SUBZONE_NO</th> <td>1</td> </tr><tr bgcolor=""> <th>SUBZONE_N</th> <td>ANAK BUKIT</td> </tr><tr bgcolor="#E3E3F3"> <th>SUBZONE_C</th> <td>BTSZ01</td> </tr><tr bgcolor=""> <th>CA_IND</th> <td>N</td> </tr><tr bgcolor="#E3E3F3"> <th>PLN_AREA_N</th> <td>BUKIT TIMAH</td> </tr><tr bgcolor=""> <th>PLN_AREA_C</th> <td>BT</td> </tr><tr bgcolor="#E3E3F3"> <th>REGION_N</th> <td>CENTRAL REGION</td> </tr><tr bgcolor=""> <th>REGION_C</th> <td>CR</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>4BAD8B2C9CEBF3F2</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20191223152313</td> </tr></table></center> 
max values  : kml_99,       <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>SUBZONE_NO</th> <td>9</td> </tr><tr bgcolor=""> <th>SUBZONE_N</th> <td>YISHUN WEST</td> </tr><tr bgcolor="#E3E3F3"> <th>SUBZONE_C</th> <td>YSSZ09</td> </tr><tr bgcolor=""> <th>CA_IND</th> <td>N</td> </tr><tr bgcolor="#E3E3F3"> <th>PLN_AREA_N</th> <td>YISHUN</td> </tr><tr bgcolor=""> <th>PLN_AREA_C</th> <td>YS</td> </tr><tr bgcolor="#E3E3F3"> <th>REGION_N</th> <td>NORTH REGION</td> </tr><tr bgcolor=""> <th>REGION_C</th> <td>NR</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>95C11920195B86C7</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20191223152313</td> </tr></table></center> 

4.5.2 Converting the Spatial* class into generic sp format

spatstat requires the analytical data in ppp object form. There is no direct way to convert a Spatial* classes into ppp object. We need to convert the Spatial classes* into Spatial object first.

The codes chunk below converts the Spatial* classes into generic sp objects.

childcare_sp <- as(childcare, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")

Let’s display the sp objects properties

childcare_sp
class       : SpatialPoints 
features    : 1925 
extent      : 11810.03, 45404.24, 25596.33, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
sg_sp
class       : SpatialPolygons 
features    : 332 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

4.5.3 Converting the generic sp format into spatstat’s ppp format

Now, we will use as.ppp() function of spatstat to convert the spatial data into spatstat’s ppp object format.

str(childcare_sp)
Formal class 'SpatialPoints' [package "sp"] with 3 slots
  ..@ coords     : num [1:1925, 1:3] 40986 28309 17829 25580 38981 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:3] "coords.x1" "coords.x2" "coords.x3"
  ..@ bbox       : num [1:3, 1:2] 11810 25596 0 45404 49301 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:3] "coords.x1" "coords.x2" "coords.x3"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +to"| __truncated__
  .. .. ..$ comment: chr "PROJCRS[\"SVY21 / Singapore TM\",\n    BASEGEOGCRS[\"SVY21\",\n        DATUM[\"SVY21\",\n            ELLIPSOID["| __truncated__
# Extract coordinates
coords <- coordinates(childcare_sp)

# Create a ppp object
childcare_ppp <- ppp(coords[, 1], coords[, 2], marks = coords[, 3])

Now, let us plot childcare_ppp and examine the different.

plot(childcare_ppp)

You can take a quick look at the summary statistics of the newly created ppp object by using the code chunk below.

summary(childcare_ppp)
Marked planar point pattern:  0 points
Average intensity 0 points per square unit
marks are numeric, of type 'double'
Summary:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
                                                

Window: rectangle = [0, 1] x [0, 1] units
Window area = 1 square unit

*** 1925 illegal points stored in attr(,"rejects") ***

4.5.4 Handling duplicated points

We can check the duplication in a ppp object by using the code chunk below.

any(duplicated(childcare_ppp))
[1] FALSE

To count the number of co-indicence point, we will use the multiplicity() function as shown in the code chunk below.

multiplicity(childcare_ppp)
integer(0)

If we want to know how many locations have more than one point event, we can use the code chunk below.

sum(multiplicity(childcare_ppp) > 1)
[1] 0

The output shows that there are 128 duplicated point events.

To view the locations of these duplicate point events, we will plot childcare data by using the code chunk below.

tmap_mode('view')
tm_shape(childcare) +
  tm_dots(alpha=0.4, 
          size=0.05)

There are three ways to overcome this problem. The easiest way is to delete the duplicates. But, that will also mean that some useful point events will be lost.

The second solution is use jittering, which will add a small perturbation to the duplicate points so that they do not occupy the exact same space.

The third solution is to make each point “unique” and then attach the duplicates of the points to the patterns as marks, as attributes of the points. Then you would need analytical techniques that take into account these marks.

The code chunk below implements the jittering approach.

childcare_ppp_jit <- rjitter(childcare_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)
any(duplicated(childcare_ppp_jit))
[1] FALSE
multiplicity(childcare_ppp)
integer(0)
sum(multiplicity(childcare_ppp) > 1)
[1] 0
tmap_mode('view')
tm_shape(childcare) +
  tm_dots(alpha=0.4, 
          size=0.05)

There are three ways to overcome this problem. The easiest way is to delete the duplicates. But, that will also mean that some useful point events will be lost.

The second solution is use jittering, which will add a small perturbation to the duplicate points so that they do not occupy the exact same space.

The third solution is to make each point “unique” and then attach the duplicates of the points to the patterns as marks, as attributes of the points. Then you would need analytical techniques that take into account these marks.

The code chunk below implements the jittering approach.

childcare_ppp_jit <- rjitter(childcare_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)
any(duplicated(childcare_ppp_jit))
[1] FALSE

4.5.5 Creating owin object

When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.

The code chunk below is used to covert sg SpatialPolygon object into owin object of spatstat.

sg_sp
class       : SpatialPolygons 
features    : 332 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

#{r} # Convert the SpatialPolygons object sg_sp to owin #sg_owin <- as(sg_sp, "owin") #

The output object can be displayed by using plot() function

#```{r}

Plot the map

#plot(sg_owin, main = “sg_owin”)

#```